Towards designing of a potential new HIV-1 protease inhibitor using QSAR study in combination with Molecular docking and Molecular dynamics simulations

Human Immunodeficiency Virus type 1 protease (HIV-1 PR) is one of the most challenging targets of antiretroviral therapy used in the treatment of AIDS-infected people. The performance of protease inhibitors (PIs) is limited by the development of protease mutations that can promote resistance to the treatment. The current study was carried out using statistics and bioinformatics tools. A series of thirty-three compounds with known enzymatic inhibitory activities against HIV-1 protease was used in this paper to build a mathematical model relating the structure to the biological activity. These compounds were designed by software; their descriptors were computed using various tools, such as Gaussian, Chem3D, ChemSketch and MarvinSketch. Computational methods generated the best model based on its statistical parameters. The model’s applicability domain (AD) was elaborated. Furthermore, one compound has been proposed as efficient against HIV-1 protease with comparable biological activity to the existing ones; this drug candidate was evaluated using ADMET properties and Lipinski’s rule. Molecular Docking performed on Wild Type, and Mutant Type HIV-1 proteases allowed the investigation of the interaction types displayed between the proteases and the ligands, Darunavir (DRV) and the new drug (ND). Molecular dynamics simulation was also used in order to investigate the complexes’ stability allowing a comparative study on the performance of both ligands (DRV & ND). Our study suggested that the new molecule showed comparable results to that of darunavir and maybe used for further experimental studies. Our study may also be used as pipeline to search and design new potential inhibitors of HIV-1 proteases.


Introduction
Human Immunodeficiency Virus (HIV) is one of the most challenging viruses in medicine, causing severe complications related to human health [1]. HIV

Chemical compounds and descriptors
Ten HIV-1 protease inhibitors have been approved by the Food and Drug Administration (FDA), but the emergence of multidrug-resistant (MDR) strains has limited long-term treatment options [19][20][21][22]; therefore, the search for new efficient drugs has become a necessity. Thirty-three new compounds were synthesized and evaluated in a previous study to determine their optimal biological activity [14] (S1 Table in S1 File).
Meanwhile, this work is based on computing various descriptors (Topological, Constitutional, Geometrical, Physicochemical, and Quantum) of the compounds mentioned above using several software packages (Gaussian, Chem3D, ChemSketch, and MarvinSketch).

Descriptive analysis
The computed descriptors must be analyzed to generate a computational model that relates the biological activity of these compounds to the structure (descriptors).
To do so, we used both methods; the first one is called Principal Component Analysis (PCA), the main purpose of which is to delete correlated descriptors, so we lower the dimension of the data representation area. The second one is a clustering method, called k-means partitioning, used to split the dataset into training sets for model generation and test sets for validation.

Statistical analysis
Multiple linear regression. Multiple linear regression analysis is a statistical technique based on several analytical independent variables called descriptors to anticipate the outcome of a response variable (biological activity); it is selected to asset a linear model relating the activity (dependent variable) to descriptors having high correlation with the response (activity) [23].
The linear model takes the form that follows: Where; Y represents the biological activity (dependent variable), a 0 is the intercept of the equation, x i is the molecular descriptors, and a i is their coefficients. Model generation. A QSAR model was generated using XLSTAT software after analyzing the data with both methods (PCA and K-means) [24], which after validation, were used to anticipate the activity of brand-new compounds that can be more efficient as HIV-1 protease inhibitors.
In order to assess the physicochemical influence of the substituents (structure/descriptors) on the biological activity, we introduced the dataset along with descriptors corresponding to the 33 compounds listed previously and their biological activities to an MLR analysis.
To choose the first-rate regression performance, we use several coefficients; r, r 2 , r 2 adj , MSE and P value [25,26], where r represents the correlation coefficient, r 2 is the coefficient of determination, r 2 adj is the coefficient adjusted for degrees of freedom, MSE is the mean squared error, and P value is the probability of Fisher statistics.
Model validation. The model generated by MLR analysis must be validated to evaluate its significance and ensure its accuracy prediction ability. In order to achieve this, we use internal and external validation.
Internal validation. Also called leave-one-out cross-validation (LOOCV), whereby one element is removed from the training set, and the remaining compounds are used to rebuild a model; then it will be returned to the training set, and another compound will be removed, the model generated will be used to predict the activity of the removed one and the cycle is repeated until all compounds have been detached one by one, in the end, a correlation coefficient Q 2 is computed [27].
External validation. Besides the internal validation, external validation is primordial; the k-means clustering method allowed us to divide the dataset into training and test sets. The second one was employed in this stage. The obtained model will be used to investigate the activities of the test set compounds, and the regression coefficient (R 2 cv ) value will be computed [28].
Applicability domain. The model was obtained based on the training set, so it is valid only with compounds with similarities as compared to those included in the training set. Therefore, new molecules must belong to the training domain. A model without an applicability domain can presume the activity of all compounds, regardless of their features, compared to those counted in the aberrant training set. So the AD is a tool to detect compounds outside the applicability domain of the obtained QSAR model and the outliers in the training set [29].

Molecular docking
Molecular Docking is an important technique used to preview the binding affinities for a vast number of small molecules, with the protease generating several conformations of the ligandprotease complex that will be ranked based on their affinity [30].
The main purpose of molecular docking study is to assess the binding energy as well as the interaction types between the ligands and the protease [31].

ADMET properties
The Absorption, Distribution, Metabolism, Elimination, and Toxicity (ADMET) properties are crucial for the effectiveness and safety of a therapeutic compound. More than 50% of practical clinical tests are unsuccessful due to the insufficiency in ADMET properties [32]. Therefore, computing ADMET properties using various servers in the drug design field can significantly shorten the probability of drug evolution failure.
These properties can be predicted using many servers, such as pkCSM [33] and Swis-sADME [34]. The obtained properties contain drug-likeness prediction based on Lipinski's rule. When compounds meet Lipinski's rule with a bioavailability score of 0.55 they will be considered as sufficiently absorbable via oral route [35,36].

Molecular dynamics simulation
Molecular dynamics simulation is the most incredible tool to predict the properties of new particles and their motion [37]. In this work, we aim to predict the dynamics information between the HIV1-protease and the proposed ligand in order to check the stability results of the docked complex [38]. For the Molecular Dynamics Simulations and MM-PBSA calculations, a similar methodology performed in a previous study was adopted [39].

Chemical compounds
A series of thirty-three compounds (inhibitors with purine base amine-acetamide as P2-ligands) synthesized and evaluated for their biological activities in previous work are the key elements in the current research; their molecular structures are listed in the ST1 Table in S1 File.

Dependent variable values
The experiment IC 50 , biological activity values, were transformed to the negative logarithm of IC 50 , using the following equation: pIC 50 = -log (IC 50 ). The results are listed in the table below ( Table 1).

Descriptors generation
Several softwares were used to compute various descriptors such as Gaussian, Chem3D, ChemSketch and MarvinSketch, but only some descriptors correlated with the activity were used in minimizing the size of the data representation space. Considering the quantum descriptors, they were investigated using DFT approach performed by Gaussian 09 program package; employing for this purpose the hybrid method B3LYP combining the Becke's threeparameter and the Lee-Yang-Parr exchange-correlation functional, using as well 6-31G (d,p) basis set, performing the optimization of the compounds geometries ultimately while all the other parameters were computed using Chem3D, Chemsketch and MarvinSketch software (S2 Table in S1 File).

Principle component analysis
Using the Principal Component Analysis, the size of the data representation space was reduced using descriptors that show a correlation coefficient with the activity higher than 0.1 in absolute value (Table 2), as well as the absence of collinearity between descriptors used to elaborate the model, was inspected by the correlation matrix.

K-Means Cluster Analysis (k-MCA)
A clustering method, called k-means partitioning, was used to cut the dataset into a training set for model generation and a test set for its validation ( Table 3). The data set is divided into five clusters. Five compounds are selected randomly, one from each cluster, to form the test set (16f, 17g, 18d, 16b and 17h), while the remaining compounds will form the training set. The last one is the key element to generate the model, and the first one was used to validate it.

Multiple linear regression (MLR)
Model generation. The model was elaborated using XLSTAT, statistical software, used as add-on for Excel. MLR equation: pIC50 = − 3-0.59*E Gap +1.27*HLC-0.033*PSA-0.015*DE Statistical parameters: R 2 = 0.66; R 2 Adj = 0.60; MSE = 0.18; P value <10 −4 ; F = 11.23 For the model above, P value is lower than 0.0001, which means that taking the risk of 0.01% by considering the null hypothesis (no effect of the descriptors on the activity) as wrong, therefore, we can assume that the model proposed includes variables with a representative amount of information. The higher values of R 2 and R 2 Adj and the lower value of MSE show that the proposed model has a higher predictive ability and reliability.
The existence of multi-collinearity among the descriptors was investigated with a parameter called variance inflation factor (VIF), the highest value is less than ten (VIF = 0.62) which further confirmed the absence of multi-collinearity problem [40,41]. The table below shows the variance inflation factor values ( Table 4).

Model interpretation:
In the proposed model, descriptors that are influencing the activity negatively are the Energy Gap (E Gap ), the Polar Surface Area (PSA) and the Dreiding Energy (DE), while only one parameter has a positive influence on the activity, which is Henry's Law Constant (HLC).
• E Gap displays a negative sign in the model, which means that increasing the activity requires minimizing E Gap value, as well as PSA and DE.
• HLC shows a positive sign in the model, allowing us to conclude that increasing the activity is achieved by increasing HLC.
To sum up, the biological activity is influenced by four variables (E Gap , HLC, PSA and DE). To increase the biological activity, E Gap must be decreased, PSA as well as the DE while HLC is increased.
Internal and external validation. The model proposed, despite its statistical parameters, must be validated following two steps:

Internal validation (Y-randomization test):
The leave-one-out cross-validation technique obtains the model's cross validation coefficient, the coefficient Q 2 LOO obtained is used as a proof of both robustness and predictive capacity of the model [42]. The given model's robustness was confirmed with a cross validation value of 0.53 (Q 2 LOO = 0.53). Y-randomization test Y-scrambling is performed on the training set; it is used to confirm that the developed model was not a result of random correlation between the biological activity and the descriptors. In this analysis, the dataset is permuted; the biological activity values were randomly distributed while the descriptors matrix was unchanged, followed by MLR analysis generating new models [43]. For each randomization and subsequent MLR analysis, we obtain a new set of values for R 2 Rand and Q 2 Rand [44] (Table ST3). If the new QSAR models have lower determination coefficient (R 2 Rand ) and leave one out determination coefficient (Q LOO 2 ) values as well for several trials (100 times in this study), we consider the proposed QSAR model as robust. Moreover, if the cRp 2 is greater than 0.5, it will be confirmed that the model is not a result of chance correlation [45,46]. For the current work, the average values of R Rand , R 2 Rand and Q 2 cv (Rand) are 0.35, 0.14 and -0.29 respectively, the cRp 2 value equals 0.60 which is higher than 0.5 (S3 Table in   random correlation between the activity and the descriptors affecting significantly the response and the developed QSAR model is robust.

External validation:
The model then must be externally validated using the test set mentioned above, in this stage, the model proposed must conclude the activities of the test set compounds in arrangement with the experimental values (Table 5), graphically presented in the figure bellow (Fig 2). The predictive ability was confirmed with a test coefficient value of 0.64 (R 2 Test = 0.64).

Applicability domain (AD)
The standardized residuals and the leverage were both jointed to illustrate the applicability domain. The Williams plot for the QSAR model is illuminated in figure below (Fig 3). The warning leverage (h*) was found to be 0.45 for the developed QSAR model. Based on the leverages, all compounds were found to be inside the defined AD.

New drugs elaboration
In order to suggest new efficient compounds, we must select from the series of compounds used in the present work, those with the highest values of pIC50 (1.37, 1.38, 1.17, 1.10) corresponding to (16a, 16f, 16j and 16k) respectively. These particles will be the object of structural modification in order to design new molecules; their descriptors' values are determined using the same tools as well as pIC 50 values predicted by MLR model proposed. Furthermore, 24 compounds candidates were designed and their parameters were computed. The leverage values (hi) were computed using Matlab software with the following equation: hi = x i T (X T X) -1 x i (i = 1, 2 . . . n) (S4 Table in S1 File).
With: x i represents the proposed compounds descriptors' matrix, X represents the test set descriptors' matrix and X T represents the transpose of the test set descriptors' matrix.
Among the 24 compounds, only one compound (16 th ) has a leverage value (h i = 0.43) lower than h* (h* = 0.45) and a biological activity higher than the known ones (pIC 50 = 1.58) (Fig 4).

ADMET properties
In the one hand, regarding Lipinski's rule, the drug-likeness of the proposed compound was verified with only one violation (MW>500) ( Table 6), which means that the proposed compound is considered as sufficiently absorbable via oral route with a bioavailability score of 0.55, in the other hand, ADMET properties predictions for the selected compound were performed using SwissADME and pkCSM web servers.
The pharmacokinetic parameters (ADMET) (absorption, distribution, metabolism, excretion, and toxicity) related to the brand-new drug are computed using pkCSM.
The absorption of the drug is primarily based on the factors that comply with; water-solubility, membrane permeability (Caco-2), intestinal absorption (human), skin permeability, p- glycoprotein. The drug distribution properties are expected from the data of volume distribution (VDss), the fraction of unbound drug, the blood-brain barrier (BBB), and central nervous system (CNS) permeability. For the biotransformation evaluation, participants of the cytochrome P450 (CYP) superfamily are selected (CYP 2D6, CYP 3A4, CYP 1A2, CYP 2C19, CYP 2C9, CYP 2D6 & CYP 3A4), while the excretion of compounds involves the total clearance of xenobiotics and renal clearance via organic cation transporter 2 (OCT 2). The toxicity of compounds is investigated using AMES toxicity; maximum tolerated dose, the human Ether-a-gogo Related Gene (hERG) potassium channel inhibition, oral rat acute toxicity, oral rat chronic toxicity, skin sensitization, T.Pyriformis toxicity and Minnow toxicity. Just a few of the important factors are mentioned in the present study, notably:

Water solubility
For the oral administrative drugs discovery, water solubility prediction is highly required. The decimal logarithm of the molar solubility in water is -3.224 (log mol/L). Considering what follows (Insoluble < -10 < poorly soluble < -6 < Moderately < -4 < soluble <

PLOS ONE
In silico designing of a potential new HIV-1 protease inhibitor -2 < very soluble < 0 < highly soluble) [47], the compound has a good solubility in water, therefore the development and the production as well of oral solid dosage is possible.

Caco2 permeability
If the predicted Papp log value is higher than 0.90 10 −6 cm/s [47], the compound is considered to have high Caco-2 permeability, for the drug candidate, it has for value 1.098 10 −6 , so we can say it has a high permeability in Caco-2.

Intestinal absorption (human)
The quantity absorbed of the drug candidate by the intestinal system is one of the major factors for oral bioavailability [48]. For the proposed compound, the intestinal absorption (human, % absorbed) seems to be 74.616%.

BBB permeability
The BBB permeability of the drug candidate has a value of -1.118 log BBB. According to the research [33], the compound is adept to cross the blood-brain barrier, if the Log BB value is higher than 0.3 and it can't cross adequately the blood-brain barrier if the log BB value is lower than -1. Therefore, the drug candidate won't be able to cross the blood-brain barrier.

CYP2D6 substrate
Drug that inhibit or compete for CYP2D6 can conduct clinical problems; this isoenzyme is highly polymorphic and is responsible for metabolizing relatively 25% of known pharmaceuticals [49]. In the current study, the drug candidate is not inhibitor of CYP2D6 enzymes.

Total Clearance
The compound has a total clearance of 0.288 log ml/min/kg, therefore, it could be excreted quickly [47].

AMES toxicity
The compound is AMES negative and test suggests that the compound could be not mutagenic [47].

hERG inhibitor
Drugs that block these HERG K+ channels are likely to cause cardiac toxicity [50]. The safe range for an ideal drug should be -5 or higher, if the value is below this level, it is predicted to cause cardiac toxicity [47].

Oral Rat Acute Toxicity (LD50)
The proposed compound is dangerous only at huge doses regarding its high LD50 value (2,259 mol/kg) [50].

Molecular docking
Molecular docking study was carried out with the aim of predicting the best conformation of the HIV-1 protease of both types (mutant and wild), on the one hand; combined to the proposed compound as a new efficient drug candidate (ND), on the other hand; combined to an FDA approved drug called Darunavir (DRV). We selected both types of the HIV-1 protease (WT and MT) as receptors. The structures of the wild type (WT) as well as the mutant type (MT) proteases were downloaded from Protein Data Bank (PDB), their PDP codes are respectively: (4LL3-Structure of wild-type HIV-1 protease in complex with Darunavir) (Fig 5) and (3TTP-Structure of multiresistant HIV-1 protease in complex with Darunavir) (Fig 6). Their original ligands were eliminated using Discovery Studio, polar hydrogens were added and the proteins were saved in PDB format, and then saved in PDBQT format using Autodock MGL Tools. The ligand proposed as a new efficient drug was earlier designed and optimized using Gaussian, then saved in PDBQT format by Autodock MGL tools (Fig 7); in addition, DRV was taken from the crystal structures downloaded from Protein Data Bank (Fig 8).
Command prompt and Vina folder were used in order to run the Docking. Different conformations of the ligand binding modes for both types were obtained with their respective binding energies (kcal/mol) after the accomplishment of the docking runs; the best pose is the one with the lowest affinity value.  The interactions between the ligands (ND & DRV) and the proteases were visualized using Discovery Studio (Table 7). Active residues interacting with the ligands (ND & DRV) are also disclosed (S5 Table in S1 File). Moreover, atoms from ligands and residues interacting with each other to form hydrogen bonds are mentioned (S8 Table in S1 File).

Molecular dynamics simulation
To evaluate the native proteins' stability (WT & MT), as well as the docked compounds' (WT-DRV, WT-ND, MT-DRV & MT-ND), a computational process is carried out through the Molecular Dynamics simulation (MD) study, allowing structural analysis at the atomic level, aiming at investigating the motion of the four complex compounds and the native proteins.
Therefore, MD simulations were administered in nine plots, with 100ns for each, using the best poses generated based on the docking study performed previously, the compounds were carried out in water simulations separately. Further, the stability analysis was performed through several techniques, namely: Root Mean Square Deviation (RMSD), Root Mean Square fluctuation (RMSF) and the Radius of Gyration (Rg).
Root Means Square Deviation (RMSD). RMSD stands for Root Means Square Deviation, it is a numerical measurement, it estimates the approximate distance between a band of atoms, mainly, backbone atoms of a protein plotted against time. The Root Means Square Deviation value is typically a measure of how much the protein's structure has been modified over time in comparison to the starting point. Further, if the RMSD of the protein presents considerable fluctuations, then no equilibrium is reached, therefore, more simulation time is required for better results.
As the RMSD plots display (Fig 17), the native proteins (WT and MT) do not show any promising stability within the simulation time especially for the wild type protease. Regarding the RMSD plots (Fig 17 (WT)) for the two complexes (WT-DRV and WT-ND), it is highly clear that these compounds are showing lower fluctuations than the native protein (WT) within the simulation time. As for the complexes (MT-ND and MT-DRV), they're showing as well lower fluctuations as compared to the native protein (MT) within the simulation time (Fig  17 (MT)). However, WT-ND and MT-ND complexes are showing promising results comparable to those of Darunavir in terms of fluctuations.
Root Means Square Fluctuation (RMSF). The Root Mean Square fluctuation (RMSF) measures the approximate deviation of a particle over time from a reference position at a specific temperature and pressure. The RMSF analysis illuminates the fluctuations of residues during the MD simulation time.
Considering the graphics, for the wild type and the mutant type proteases for both chains (Fig 18)  Darunavir are significantly influencing the fluctuations of the proteins' residues in most regions.

Radius of gyration (Rg).
The radius of gyration is an interesting parameter as well to investigate the motion of a protein as well as its stability; it describes the compactness of the protein during the simulation time.
For the Mutant Type protease (Fig 19 (MT)), the radius of gyration of the complex compound MT-ND is higher in value as compared to the MT native protein and the complex compound in presence of DRV, causing eventually higher flexibility of the compound MT-ND. For the Wild Type protease (Fig 19 (WT)), the plots show that the complex compound WT-ND reveals more compactness with lower radius of gyration values as compared to the complex compound WT-DRV and the WT native protein within the simulation time, inducing less flexibility, which means higher potential of stability for the complex WT-ND.
Hydrogen bonds. Hydrogen bonds are primordial in drug specificity and stability, so the determination of H-bond number in complex compounds is essential to check its contribution to the overall stability of each system and further conduct a comparative study including all complex compounds in question.
The figure (Fig 20) shows that during the MD simulation period (100ns), the complex MT-ND's graph is showing up to seven hydrogen bonds by the end of the simulation time, while the MT-DRV complex compound's graph is showing a few hydrogen bonds during the first 40ns as compared to MT-ND, then significantly increasing at 60ns displaying ten hydrogen bonds then decreasing to seven by the end of the simulation time (Fig 20 (MT)). In contrast, for the complex compound WT-DRV, the number of hydrogen bonds is consistently decreasing from 8 to 5 while the compound WT-ND displays up to five hydrogen bonds with no significant decrease compared to the WT-DRV compound during the simulation time (Fig  20 (WT)). We can conclude that whether the wild type or the mutant type proteases, when docked to the new drug, the number of hydrogen bonds is likely to be the same with no significant change as compared to the complex compounds with Darunavir that shows a decreasing number of hydrogen bonds during the simulation time.
Hydrophobic interactions. Hydrophobic interactions are non-bonded interactions between the protein and the ligand, which play a major role in the stability of complexes.
As shown below, considering the wild type protease (Fig 21 (WT), both complexes WT-DRV and WT-ND show highly similar numbers of hydrophobic interactions during the simulation time. In contrast, for the mutant type protease (Fig 21 (MT), the complexes MT-DRV and MT-ND, the number of hydrophobic interactions for the complex compound MT-DRV is significantly higher than the number of hydrophobic interactions for the complex compound with the new drug MT-ND. We can conclude that for the wild-type protease, the new drug significantly competes with Darunavir, displaying similar numbers of hydrophobic interactions at every 20 ns of the simulation time. However, Darunavir is showing highly promising results for the mutant-type protease compared to the new drug in terms of hydrophobic interactions. Solvent Accessible Surface Area (SASA). The accessible surface area (ASA) or solventaccessible surface area (SASA) is the surface area of a biomolecule that is accessible to a solvent.
Based on the graphics (Fig 22), the new drug, when combined to the wild type protease, is showing promising results regarding the significant decrease of the ASA values since 40ns to the end of the simulation time (Fig 22 (WT)), but for the mutant type, the ASA values are not promising on the ground that the graphics are displaying increasing values starting from 60ns of the simulation time (Fig 22 (MT)).
We can conclude that the new drug is comparable to Darunavir during the last 30ns of the simulation time for the wild type protease while no possible competition is investigated for the mutant type on the ground that the graphic is showing significant ASA values for the complex MT-ND as compared the MT-DRV mainly during the last 40ns of the simulation time.
Binding free energy calculation. Molecular dynamics simulations were used to calculate binding free energy using the MM-PBSA method. Snapshots were extracted at every 1 ns of stable intervals from 70-100 ns MD trajectory. The binding free energy and its corresponding component obtained from the MM-PBSA calculations are listed ( Table 8).
The results indicate that for both wild and mutant type protease, Darunavir is showing a binding affinity of -173.323 kJ/mol and -190.868 kJ/mol, respectively, which is slightly higher than the New Drug (-170.903 kJ/mol and -187.521 kJ/mol, respectively). van der Waals, Electrostatic and SASA energy played a crucial role in binding energy and complex stability. In contrast, polar solvation energy has an opposite effect causing binding energy to depend on its unfavorable positive value. Among different energy terms, the contribution of van der Waals energy towards total binding energy is superior.
Compilation of the data demonstrated that although the binding of Darunavir to both wild and mutant HIV protease is better, the binding of the new drug is comparable to that of Darunavir in both wild and wild mutant type. This is illustrated by the different analyses that have been used so far. Thus, the new drug may also be considered a potential inhibitor against multi-drug resistant HIV and may be tested experimentally.

Conclusion
Various softwares have been used in this study in order to generate a reliable model relating the biological activity of new HIV-1 protease inhibitors to their physicochemical parameters. The generated model showed a high predictability efficiency regarding its statistical parameters. The applicability domain was also generated to frame the workspace (only compounds with features with greater similarity to those included in the training set can be used). Regarding the proposed model, the biological activity of the new HIV-1 protease inhibitors can be increased by increasing the three variables' values; the Energy Gap (E Gap ); the Polar Surface Area (PSA) and the Dreiding Energy (DE) (positively related to the activity), and decreasing the Henry's Law Constant value (negatively related to the activity). A new drug was proposed based on the model generated with a biological activity higher than the known drug compounds' activities. Afterwards, the molecular docking study was performed on the wild-type and the mutant-type HIV-1 proteases to predict the best conformation displayed by two ligands, the New Drug and Darunavir as an approved FDA drug. Moreover, molecular dynamics simulation was performed to study the stability of the complexes (WT-DRV, WT-ND, MT-DRV & MT-ND); results disclosed some interesting results related to the new drug, therefore, the new drug may be considered as a potential inhibitor against multi-drug-resistant (MDR) strains of HIV-1 protease (PR) and may be tested experimentally.